home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_8 / quadmin.m < prev   
Encoding:
Text File  |  1994-06-05  |  2.3 KB  |  99 lines  |  [MATF/MATL]
  1. function [p,yp,dp,dy,P] = quadmin(f,a,b,delta,epsilon)
  2. % [p,yp,dp,dy,P] = quadmin(f,a,b,delta,epsilon)
  3. % [p,yp,dp,dy] = quadmin(f,a,b,delta,epsilon)
  4. % Search for a minimum, using quadratic interpolation.
  5. % f  is the function, input.
  6. % a  is the left  endpoint, input.
  7. % b  is the right endpoint, input.
  8. % delta   is the tolerance for the abscissas, input.
  9. % epsilon is the tolerance for the ordinates, input.
  10. % p  is the abscissa of the minimum, output.
  11. % yp is the ordinate of the minimum, output.
  12. % dp is the error bound for  p, output.
  13. % dy is the error bound for yp, output.
  14. % P  is the vector of iterations, output.
  15. p0 = a;
  16. maxj = 20;
  17. maxk = 30;
  18. big = 1e6;
  19. err = 1;
  20. k = 1;
  21. P(k) = p0;
  22. cond = 0;
  23. h = 1;
  24. if (abs(p0)>1e4), h = abs(p0)/1e4; end
  25. while (k<maxk & err>epsilon & cond~=5)
  26.   f1 = (feval(f,p0+0.00001)-feval(f,p0-0.00001))/0.00002;
  27.   if (f1>0), h = -abs(h); end
  28.   p1 = p0 + h;
  29.   p2 = p0 + 2*h;
  30.   pmin = p0;
  31.   y0 = feval(f,p0);
  32.   y1 = feval(f,p1);
  33.   y2 = feval(f,p2);
  34.   ymin = y0;
  35.   cond = 0;
  36.   j = 0;
  37.   while (j<maxj & abs(h)>delta & cond==0)
  38.     if (y0<=y1),
  39.       p2 = p1;
  40.       y2 = y1;
  41.       h = h/2;
  42.       p1 = p0 + h;
  43.       y1 = feval(f,p1);
  44.     else
  45.       if (y2<y1),
  46.         p1 = p2;
  47.         y1 = y2;
  48.         h = 2*h;
  49.         p2 = p0 + 2*h;
  50.         y2 = feval(f,p2);
  51.       else
  52.         cond = -1;
  53.       end
  54.     end
  55.     j = j+1;
  56.     if (abs(h)>big | abs(p0)>big), cond=5; end
  57.   end
  58.   if (cond==5),
  59.     pmin = p1;
  60.     ymin = feval(f,p1);
  61.   else
  62.     d = 4*y1-2*y0-2*y2;     %Start of a long block:
  63.     if (d<0),
  64.       hmin = h*(4*y1-3*y0-y2)/d;
  65.     else
  66.       hmin = h/3;
  67.       cond = 4;
  68.     end
  69.     pmin = p0 + hmin;
  70.     ymin = feval(f,pmin);
  71.     h = abs(h);
  72.     h0 = abs(hmin);
  73.     h1 = abs(hmin-h);
  74.     h2 = abs(hmin-2*h);
  75.     if (h0<h),  h = h0;   end
  76.     if (h1<h),  h = h1;   end
  77.     if (h2<h),  h = h2;   end
  78.     if (h==0),  h = hmin; end
  79.     if (h<delta), cond=1; end
  80.     if (abs(h)>big | abs(pmin)>big), cond=5; end
  81.     e0 = abs(y0-ymin);
  82.     e1 = abs(y1-ymin);
  83.     e2 = abs(y2-ymin);
  84.     if (e0~=0 & e0<err), err = e0; end
  85.     if (e1~=0 & e1<err), err = e1; end
  86.     if (e2~=0 & e2<err), err = e2; end
  87.     if (e0~=0 & e1==0 & e2==0), error=0; end
  88.     if (err<epsilon), cond=2; end
  89.     p0 = pmin;
  90.     k = k+1;
  91.     P(k) = p0;
  92.   end                           % End of the long block.
  93.   if (cond==2 & h<delta), cond=3; end
  94. end
  95. p = p0;
  96. dp = h;
  97. yp = feval(f,p);
  98. dy = err;
  99.